{
 "metadata": {},
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Spherical Bessel Zeros\n",
      "======================\n",
      "\n",
      "It may be useful to find out the zeros of the spherical Bessel\n",
      "functions, for instance, if you want to compute the eigenfrequencies of\n",
      "a spherical electromagnetic cavity (in this case, you'll need also the\n",
      "zeros of the derivative of (r\\*Jn(r))).\n",
      "\n",
      "The problem is that you have to work out the ranges where you are\n",
      "supposed to find the zeros.\n",
      "\n",
      "Happily, the range of a given zero of the n'th spherical Bessel\n",
      "functions can be computed from the zeros of the (n-1)'th spherical\n",
      "Bessel function.\n",
      "\n",
      "Thus, the approach proposed here is recursive, knowing that the\n",
      "spherical Bessel function of order 0 is equal to sin(r)/r, whose zeros\n",
      "are well known.\n",
      "\n",
      "This approach is obviously not efficient at all, but it works ;-).\n",
      "\n",
      "A sample example is shown, for the 10 first zeros of the spherical\n",
      "Bessel function of order 5 (and the derivative of (r\\*J5(r))), using\n",
      "[:Cookbook/Matplotlib: matplotlib]."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#! /usr/bin/env python\n",
      "\n",
      "### recursive method: computes zeros ranges of Jn(r,n) from zeros of Jn(r,n-1)\n",
      "### (also for zeros of (rJn(r,n))')\n",
      "### pros : you are certain to find the right zeros values;\n",
      "### cons : all zeros of the n-1 previous Jn have to be computed;\n",
      "### note : Jn(r,0) = sin(r)/r\n",
      "\n",
      "from scipy import arange, pi, sqrt, zeros\n",
      "from scipy.special import jv, jvp\n",
      "from scipy.optimize import brentq\n",
      "from sys import argv\n",
      "from pylab import *\n",
      "\n",
      "def Jn(r,n):\n",
      "  return (sqrt(pi/(2*r))*jv(n+0.5,r))\n",
      "def Jn_zeros(n,nt):\n",
      "  zerosj = zeros((n+1, nt), dtype=Float32)\n",
      "  zerosj[0] = arange(1,nt+1)*pi\n",
      "  points = arange(1,nt+n+1)*pi\n",
      "  racines = zeros(nt+n, dtype=Float32)\n",
      "  for i in range(1,n+1):\n",
      "    for j in range(nt+n-i):\n",
      "      foo = brentq(Jn, points[j], points[j+1], (i,))\n",
      "      racines[j] = foo\n",
      "    points = racines\n",
      "    zerosj[i][:nt] = racines[:nt]\n",
      "  return (zerosj)\n",
      "\n",
      "def rJnp(r,n):\n",
      "  return (0.5*sqrt(pi/(2*r))*jv(n+0.5,r) + sqrt(pi*r/2)*jvp(n+0.5,r))\n",
      "def rJnp_zeros(n,nt):\n",
      "  zerosj = zeros((n+1, nt), dtype=Float32)\n",
      "  zerosj[0] = (2.*arange(1,nt+1)-1)*pi/2\n",
      "  points = (2.*arange(1,nt+n+1)-1)*pi/2\n",
      "  racines = zeros(nt+n, dtype=Float32)\n",
      "  for i in range(1,n+1):\n",
      "    for j in range(nt+n-i):\n",
      "      foo = brentq(rJnp, points[j], points[j+1], (i,))\n",
      "      racines[j] = foo\n",
      "    points = racines\n",
      "    zerosj[i][:nt] = racines[:nt]\n",
      "  return (zerosj)\n",
      "\n",
      "n = int(argv[1])  # n'th spherical bessel function\n",
      "nt = int(argv[2]) # number of zeros to be computed\n",
      "\n",
      "dr = 0.01\n",
      "eps = dr/1000\n",
      "\n",
      "jnz = Jn_zeros(n,nt)[n]\n",
      "r1 = arange(eps,jnz[len(jnz)-1],dr)\n",
      "jnzp = rJnp_zeros(n,nt)[n]\n",
      "r2 = arange(eps,jnzp[len(jnzp)-1],dr)\n",
      "\n",
      "grid(True)\n",
      "plot(r1,Jn(r1,n),'b', r2,rJnp(r2,n),'r')\n",
      "title((str(nt)+' first zeros'))\n",
      "legend((r'$j_{'+str(n)+'}(r)$', r'$(rj_{'+str(n)+'}(r))\\'$'))\n",
      "plot(jnz,zeros(len(jnz)),'bo', jnzp,zeros(len(jnzp)),'rd')\n",
      "gca().xaxis.set_minor_locator(MultipleLocator(1))\n",
      "# gca().xaxis.set_minor_formatter(FormatStrFormatter('%d'))\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bessph_zeros_rec 5 10"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}